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We present a generalization of dissipative particle dynamics that includes shear forces between 
particles. The new algorithm has the same structure as the (isothermal) smoothed particle dy- 
namics algorithm, except that it conserves angular momentum and includes thermal fluctuations 
consistently with the principles of equilibrium statistical mechanics. This clarifles the connection of 
dissipative particle dynamics with numerical resolution algorithms of the macroscopic Navier-Stokes 
equations. 

Dissipative particle dynamics (DPD) is an ofF-lattice simulation technique that has been introduced by Hooger- 
brugge and Koelman in order to address hydrodynamic problems in complex fluids [0J|]. The technique has received 
substantial theoretical support and it has been successfully applied in a large variety of systems including flow 
in porous media 1^ , colloidal suspension |^,|| , dilute polymeric suspensions , and immiscible binary fluids |^ . 

The idea behind DPD is to simulate a Newtonian fluid (as, for example, the Newtonian solvent in colloidal or 
polymeric suspensions) in terms of mesoscopic "lumps" or "droplets" of fluid, named dissipative particles It 
is postulated that these dissipative particles interact with each other with a pair-wise conservative potential, with 
dissipative forces that depend on the relative approaching velocity of the particles, and with random forces that 
satisfy a fluctuation-dissipation theorem |^ . Newton's third law is satisfied and the total momentum of the system is 
conserved, although energy is not |4|. This implies that there are local conservation equations for mass and momentum, 
i.e. the system behaves hydrodynamically at long times and distances ^]. The advantage of DPD over conventional 
molecular dynamics relies on the fact that DPD is a coarse-grained technique which captures the gross features of 
mesoscopic portions of fluid. The microscopic details, which are computationally expensive and not even interesting, 
are averaged out in DPD. 

The spirit of DPD turns out to be quite similar to that of smoothed particle dynamics SPD . This technique 
was originally intended to simulate astrophysical non- viscous flows and has been recently been applied in a variety of 
non- viscous , [0 and viscous problems ||l2| , . SPD consists on a discretization of the Navier-Stokes equations 
in a Lagrangian moving grid with the aid of a weight function. In this way, the nodes of the grid can be identified 
as "smooth particles" interacting through prescribed laws of force, thus allowing to solve Navier-Stokes equations 
with molecular dynamics codes. The Lagrangian nature of the technique makes it very appropriate to study flows in 
complex geometries (like those appearing in colloidal suspensions) because there is no need of costly recalculations 
of the mesh as the boundary conditions evolve. Actually, the very dynamics takes care of it. Unfortunately, there 
is at present no implementation of the thermal fluctuations present in fluids at mesoscopic scales and which are the 
responsible for the Brownian motion of small suspended objects. It is not clear that the fluctuations that appear as 
a consequence of the discrete nature of the technique are compatible with the principles of statistical mechanics. In 
particular, that they obey a fluctuation-dissipation theorem. In other words, there is yet no implementation of SPD 
for fluctuating hydrodynamics p[ . 

In this letter we present a model of fluid particles that interact with dissipative forces that, besides the dissipative 
force of the original DPD algorithm, include shear forces between fluid particles. We regard this as a more realistic 
model in view of the fact that the proposed algorithm coincides in structure with the SPD algorithm. Therefore, it is 
expected that an even more reduced number of particles already reproduce the hydrodynamic behavior of the system. 
We also formulate the random forces between fluid particles in such a way that the distribution function of velocities 
is Gaussian, as predicted by equilibrium statistical mechanics. In this way, this work represents a generalization of 
SPD that includes thermal fluctuations. This opens up the possibility of applying a technique closely related to SPD 
to the study of complex fluids. 

A shear force between particles i,j is proportional to the relative velocities Vjj of both particles, whereas the 
dissipative force in the original DPD algorithm is proportional to the approaching velocity (e^ •Vij)ey , where e^j is 
the unit vector in the joining line of both particles. The initial motivation for modifying the original algorithm of 
DPD by introducing shear forces was the identification of an elementary motion between dissipative particles that 
produces no force in that algorithm. If a dissipative particle is orbiting in a circumference around a reference particle, 
it will not exert any force on this particle. Nevertheless, on simple physical grounds one expect that the motion of 



*e-mail: pep@flsfun.uned.es 



2 



the dissipative particle must drag in some way the reference particle. This is taken into account through the shear 
forces in the fluid particle model presented in this letter. We note, however that this relative motion might produce 
a drag even in the original DPD algorithm if many DPD particles are involved simultaneously. The same is true for 
a purely conservative molecular dynamics simulation. The point is, of course, that the effect is already captured with 
a much smaller number of particles in the fluid particle model. 

The shear forces are not central and therefore angular momentum is not conserved. We restore angular momentum 
conservation by including in the description a spin variable. If one thinks of the fluid particles as mesoscopic portions 
of fluid, this spin variable has a sounded physical interpretation: it describes the angular momentum of the molecules 
that constitute the fluid particle with respect to the center of mass of the fluid particle. 



I. THE FLUID PARTICLE MODEL 



The fluid particle model is defined by N identical particles of mass m and moment of inertia /. The state of the 
system is characterized by the positions r^, velocities v^, and angular velocities uJi of each particle. We do not include 
here an internal energy variable and the resulting algorithm, like DPD, will not conserve energy locally. This may be 
a minor problem when one is interested only in rheological properties. 

The equations of motion of the system are given by 



= V 

V, 



m ^ 

where Fy , are the force and torque that particle j exerts on particle i. We require that the forces satisfy Newton's 
third law, Fy = — F^i, in such a way that the total linear momentum P = mvi is a dynamical invariant, P = 0. 
In addition, we assume that the torques in (]l|) are given by Ny = — x Fy /2 and one checks immediately that the 
total angular momentum J = x Pi + i^i) is conserved exactly, J = 0. 

We model the force F^ between fluid particles according to F^j = F^ + F?J + F^ + Fy where the different 
contributions are given by 



F^. = -7mM^(r,y)-v, 
F« = -7mM«(r,,).( 



X \uJi + 



F,,dt = am (^A(ry )dW;- + S(ry )-tr[rfWy]l + (7(r„)dW^ j -e,, (2) 

The first contribution Ffj is a repulsive conservative force derived from a soft potential V{r). If only this force is 
present we have the version of SPD for non-viscous flows studied extensively in Rcfs. {i.e. a MD simulation). 
The second contribution F^ is a friction force that depends on the relative translational velocities Vy = — . The 
dimensionless matrix M-^(ry) is the most general matrix that can be constructed out of the vector = — r^, 
this is M^(ry) = A(ry)l + B{r,j) eijGij where 1 is the unit matrix, e^j — Yij/vij is the unit vector joining the 
particles, = \vij \ and the functions A{r) and B{r) provide the range of the force. The friction coefficient 7 has 
been introduced as an overall factor for convenience and has dimensions of inverse of time. Fj]- is the sum of a shear 
force —^mA{r)\-ij and a central force —jmB{r){eij ■Vij)eij . The rotational contribution F^ in (||) is given also in 
terms of a dimensionless matrix M-'^ ~ M^(rij) = C(ry )1 + D{rij)eijeij where C(r), D{r) are scalar functions (even 
though the D{r) contribution to F^ is zero we maintain this term to keep the analogy between and M-^). Note 
that if particles i and j were spheres of radius rij/2 in contact and spinning with angular velocities oJi^LOj the relative 
velocity at the "surface" of the spheres would be equal to x (w^ + Wj). Then F-'^ gives a friction force between the 
spheres proportional to this relative velocity. This force produces an "engaging" effect in which neighbour particles 
prefer to spin in opposite senses. 

The last contribution in (^) is a velocity-independent stochastic force which is inspired by the tensorial structure of 
the random forces that appear in the fluctuating hydrodynamics theory [n4[ . cr is a parameter governing the overall 
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noise amplitude, the scalar functions A{r), B{r), C{r) define the range of the random force, and we have introduced 
the following symmetric, antisymmetric and traceless symmetric random matrices 



dW: 



dW. 



^3 

Afiv 



dW 



dW: 



dW^^] 



dW^^] 



tr[dWf,]l 



(3) 



Here, D is the physical dimension of space and dW^-" is a matrix of independent Wienner increments which is assumed 



to be symmetric under particle interchange dW 



This symmetry will ensure momentum conservation be- 



cause Fij = —Fji. The matrix dW^^ is an infinitesimal of order 1/2, and this is summarized in the Ito mnemotechnical 

rule dW^jf dWJjv' = [5ij5i'j' + 5iji5ji']5^^u5^j,'u'dt. From this stochastic property, one derives straightforwardly the 
following rules from the different parts 

tr[dWii,]tr[dWjjv] = [5ij8,,j, +5,f5j,,]Ddt 
tr[dWi,.]dWf-, = tr[dW,, 



dwl'^^'dw'^;,"' 



1a 



dt 



]dWfj, = dw^r'dwff' 







(4) 



These expressions show that the traceless symmetric, the trace and the antisymmetric parts are independent stochastic 
processes. The apparently complex structure of the random force is required in order to be consistent with the tensor 
structure of the dissipative friction forces. This will become apparent when considering the associated Fokker-Planck 
equation and requiring that it has a proper equilibrium ensemble. 

The force F.y is the most general force that can be constructed out of the vectors ri,rj,'Vi,-Vj,uji,u!j and satisfies 
that : 1) it is invariant under translational and Galilean transformations and transforms as a vector under rotations; 
2) it is linear in the linear and angular velocities. This linearity is required in order to be consistent with the Gaussian 



distribution of velocities at equilibrium, as we will show later; 3) it satisfies Newton's third law Fij 
therefore, the total linear momentum is a conserved quantity of the system. 



-Fji and. 



II. FOKKER-PLANCK EQUATION AND EQUILIBRIUM STATE 



The equations of motion (|l|) are Langevin equations which have associated a mathematically equivalent Fokker- 
Planck equation The FPE governs the distribution function p{r, v, lu] t) that gives the probability density that the 
N particles of the system have specified values for the positions, velocities and angular velocities. Following standard 
pa although quite tedious procedures, the FPE is given by 



dtp{r, V, uj- 1) = [L^ + + L^-] p{r, v, lu; t) 



(5) 



where 
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Here, the matrix Ty is given by 
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The steady state solution of equation (||), dtp = 0, gives the (unique) equihbrium distribution p'^*. We now consider 
the conditions under which the steady state solution is the Gibbs canonical ensemble 



p'"^{r,v,u}) 



^} + V[r) /ksT} 



(8) 



where V is the potential function that gives rise to the conservative forces F*-^, fc^ is Boltzmann's constant, T is the 
equilibrium temperature and Z is the normalizing partition function. We note that the velocity and angular velocity 
hydrodynamic fields are Gaussian variables at equilibrium and, therefore, one expects that the distribution function 
of the discrete values of theses fields is also Gaussian. 

The canonical ensemble is the equilibrium solution for the conservative system, i.e. L'-^ p^'' = 0. If, in addition, 
the following equations are satisfied ^JjP^'^ = ^fjP'^'^ — then we will have Lp^'' — and the Gibbs equilibrium 
ensemble will be the unique stationary solution of the dynamics. These equations will be satisfied if the detailed 
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is satisfied and also M^(rij) = M-^(r,y) 



This implies 
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(9) 



We observe, therefore, that the initial hypothesis for the tensorial structure of the dissipative and random forces was 
correct and consistent with equilibrium statistical mechanics. 

The structure of the dissipative forces postulated in the fluid particle model (disregarding angular variables) is es- 
sentially the same as the structure of the viscous forces in the (isothermal) SPD algorithm. By using the discretization 
of Takeda et al. the correspondence is 
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(C+l) \w"(r) 



W'{r) 



W'{r) 



(10) 



where poi '^0 are the equilibrium pressure and number density, and W{r) is the bell-shaped weight function used in the 
discretization of the Navier-Stokes equation (the assumption that the density of all particles is almost constant has 
been taken). The SPD algorithm does not conserve angular momentum and does not incorporate thermal fluctuations. 
The first issue can be reduced at the cost of increasing the resolution (i.e. by increasing the computational cost of 
the simulations jl^). Regarding the second issue, we have formulated in this letter how to introduce consistently the 
thermal noise (i.e. by selecting A{r), B{r), C{r) satisfying (^). However, we note a serious problem in the expression 
for the scalar function A{r) in terms of the weight function W{r) in ([l0|). The left hand side of the second equation 
in (10) is negative for some values of r for a bell-shaped weight function W{r). This is unacceptable in view of Eqn. 
(^. The smoothed particle dynamics algorithm does not allow, then, for a consistent introduction of thermal noise, 
at least in the form presented in Rcf. [|l2| . 

In summary, we have proposed a fluid particle model by introducing shear forces and spin into the original DPD 
algorithm. The model has the correct equilibrium state and has a structure similar to SPD, but with angular 
momentum conservation and correct thermal fluctuations. In this way, the connection between DPD and SPD is 
clarified. 

This work has been partially supported by a DGICYT Project No PB94-0382 and by E.G. contract ERB-GHRXGT- 
940546. 
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